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1.1 Introduction 

The properties of systems near quantum critical points (QCPs) have been studied 
extensively in recent years |[Tl[2l. A QCP is a point across which the symmetry of 
the ground state of a quantum system changes in a fundamental way; such a point 
can be accessed by changing some parameter, say A, in the Hamiltonian governing 
the system. The change in the ground state across a QCP is mediated by quantum 
fluctuations. Unlike conventional thermal critical points, thermal fluctuations do not 
play a crucial role in such transitions. Similar to its thermal counterparts, the low 
energy physics near a QCP is associated with a number of critical exponents which 
characterize the universality class of such a transition. Amongst these exponents, the 
dynamical critical exponent z provides the signature of the relative scaling of space 
and time at the transition and has no counterpart in thermal phase transitions. The 
other exponent which is going to be important for the purpose of this review is the 
well-known correlation length exponent v. These exponents are formally defined as 
follows. As we approach the critical point at A = , the correlation length diverges 
as ^ ~ |A — Ac|^^, while the gap between the ground state and first excited state van- 
ishes as AE ~ ^ ^ |A — A^l^'^. Exactly at the critical point A = A^, the energy of the 
low-lying excitations vanishes at some wave number ko as © ~ |k — ko|''. The crit- 
ical exponents are independent of the details of the microscopic Hamiltonian; they 
depend only on a few parameters such as the dimensionality of the system and the 
symmetry of the order parameter. These features render the low energy equilibrium 
physics of a quantum system near a QCP truly universal. 

In contrast to this well-understood universality of the equilibrium properties 
of a system near a QCP, relatively few universal features are known in the non- 
equilibrium behavior of a quantum system. Initial studies in this field aimed at un- 
derstanding the near-equihbrium finite temperature dynamics near a quantum critical 
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point using the Boltzman equation approach |[3]- Such a dynamics is useful in making 
contact with experiments which are always carried out at finite temperature. More- 
over, the excitations near a quantum critical critical point with a non-zero value of rj 
do not have a simple pole structure like that of the conventional quasiparticle excita- 
tions of condensed matter systems; this property makes such a dynamics interesting 
in its own right. 

More recently, significant theoretical llH |5] |6] |7] |8] and experimental [|9] endeav- 
ors have focussed on out-of-equilibrium dynamics of closed quantum critical sys- 
tems. On the experimental front, it has been possible, in ultracold atom systems, to 
gain unprecedented control over the measurement of out-of-equilibrium properties 
of quantum systems [9]. On the theoretical front, such studies can be broadly classi- 
fied into two distinct categories. The first type involves a study of the time evolution 
of a quantum system after a rapid quench through a quantum critical point. Such 
a study yields information about the order parameter dynamics across a quantum 
critical point. It turns out that such a dynamics exhibits a universal signature of the 
quantum critical point crossed during the quench. The second type involves a study 
of defect production during slow non-adiabatic dynamics through a quantum crit- 
ical point. Such a defect production mechanism was first pointed out for dynamics 
through thermal critical points in Refs. ifTolfTTI . For a slow enough quenches through 
quantum critical points, the density of defects produced are known to depend on z 
and V which characterize the critical point lfT2l [131 [141 . 

Quantum communication in spin systems has also been a subject of intense study 
recently. Following the seminal work in Ref. {Tsj . a tremendous amount of theoret- 
ical effort has been put in to understand the nature of qubit transfer through one- 
or multi-dimensional spin systems 1T61I . One of the major goals of such studies is 
to characterize the fidelity of the transfer of a qubit across such a spin system. The 
maximization of both the fidelity and the speed of transfer, in moving a qubit through 
a spin chain, is an issue of great interest in such studies. 

In this article, we will review some studies of sudden and slow zero temperature 
non-equilibrium dynamics of closed quantum systems across critical points. In Sect. 
11.21 we consider a sudden quench across a quantum critical point. We study the 
order parameter dynamics of one-dimensional ultracold atoms in an optical lattice 
in Sect. ll.2.Tl and of the infinite range ferromagnetic Ising model in Sect. 11.2.21 We 
demonstrate that the dynamics shows universal signatures of the QCP across which 
the system is quenched. In Sect. 11.31 we discuss defect production for slow non- 
adiabatic dynamics; typically, we find that the density of defects scales as an inverse 
power of the quench time T, where the power depends on the dimensionality d of the 
system, and the exponents z and v. In Sect. II. 3. 11 we discuss the time evolution of 
the system across a quantum critical surface; we find that the defect scaling exponent 
in this case depends on the dimensionality of the critical surface. This is confirmed 
by a study of defect production in the Kitaev model, which is an exactly solvable 
model of spin-l/2's on a honeycomb lattice. In Sect. 11.3.21 we study the effect of 
quenching across a QCP in a non-linear way; we find that the defect scaling exponent 
also depends on the degree of non-linearity. We illustrate these ideas by studying 
two exactly solvable spin- 1/2 models in one dimension. In Sect. 11.3.31 we discuss 
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a number of experimental systems where our results on defect scaling can possibly 
be checked. Finally, in Sect. |1.4[ we show that non-equilibrium dynamics, in one- 
and two-dimensional Heisenberg spin models, can be engineered to maximize the 
fideUty and speed of the transfer of qubits. 



1.2 Quench dynamics 

1.2.1 Ultracold atoms in an optical lattice 

In this section, we shall study a system of ultracold spinless bosons in a one- 
dimensional (ID) optical lattice in the presence of a harmonic trap potential fVl]. 
We will restrict ourselves to the Mott phase of the bosons and will study their re- 
sponse to a shift in the position of the trap potential. Such a shift acts as an effective 
"electric field" for the bosons whose Hamiltonian is given by Iil7i 

^ - -tL{bJbj + b]b,') + I X^«K«<-1) - E £eT,«,-, (1.1) 

ij i i 

where ij represents pairs of nearest neighbor sites of the optical lattice, n, = b^bi is 
the number operator for the bosons, r, are the dimensionless spatial coordinates of 
the lattice sites (the lattice spacing is unity), e is a unit vector in the direction of the 
applied electric field, and the effective electric field E (in units of energy) can be 
deduced from the shift a of the center of the trap as E = —adVtrap{x)/dx. In what 
follows, we will restrict ourselves to \U — E\,t <^ E,U. We note that such a regime 
has been achieved in experiments |9|. 

In the presence of such an electric field, our classical intuition suggests that all the 
bosons would gather in the last site of the ID chain thereby minimizing their energy. 
However, this does not happen for two reasons. First, the bosons are interacting and 
a state where all the bosons are in a single site leads to a huge interaction energy 
cost. But more importantly, even non-interacting bosons (or in the parameter regime 
£■ 3> t/«o for interacting bosons) do not exhibit this behavior. To understand this, 
we note that when U — 0, is simply the Wannier-Stark Hamiltonian whose wave 
functions, in the limit of strong electric fields (f <C E), are well localized Bessel 
functions. Thus for E ^ t, the bosons remain localized in their respective lattices. It 
turns out that for realistic optical lattices where interband energy spacings are large 
compared to both U and E, the Zener tunneling time, i.e., the time taken by the 
bosons to reach the final ground from this metastable Mott state, is of the order of 
milliseconds and is larger than the system lifetime [9J- Our strategy will therefore be 
to start from the parent Mott state of these localized bosons, identify the complete 
set of states resonantly coupled to this parent state, obtain the effective Hamiltonian 
within the subspace of these states, and determine its spectrum and correlations. This 
effective Hamiltonian is expected to describe the low energy behavior of the system. 

The parent Mott state and its resonant dipole excitations are shown in Figs. 11.11 
and ll.2l L17J . A dipole here consists of a bound pair of hole at site / and an additional 
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Fig. 1.1. Schematic representation of the parent Mott insulating state with no = 2. Each well 
represents a local minimum of the optical lattice potential — we number these as 1-5 from 
the left. The potential gradient leads to a uniform decrease in the on-site energy of an atom 
as we move to the right. The grey circles are the di bosons of Eq. il.2\ . The vertical direction 
represents increasing energy: the repulsive interaction energy between the atoms is realized 
by placing atoms vertically within each well, so that each atom displaces the remaining atoms 
upwards along the energy axis. We have chosen the diameter of the atoms to equal the po- 
tential energy drop between neighboring wells — this corresponds to the condition U = E. 
Consequently, a resonant transition is one in which the top atom in a well moves horizontally 
to the top of a nearest neighbor well; motions either upwards or downwards are non-resonant. 

particle at its neighboring / + 1 site. We note that the dipole excitations cost an energy 
U — E and hence become energetically favorable when the electric field exceeds 
the interaction energy. However, once a dipole forms between two adjacent sites, 
these sites cannot participate in the formation of another dipole since the resultant 
state lies out of the resonant subspace |17|. This leads to a constraint on the dipole 
number on any given link I connecting two sites, namely, < 1 . Similar reasoning, 
elaborated in Ref. ifTTl , shows that there can be at most one dipole on two adjacent 
links: w^n^+i = 0. The effective Hamiltonian of these dipoles can be written in terms 
of the dipole annihilation and creation operators di and d\ as 

= -f\/«o(«o+l) + + {U-E)Y,dldt. (1.2) 

Note that the presence of boson hopping leads to non-conservation of the dipole 
number since such a hopping can spontaneously create or destroy dipoles on a given 
link. Also, needs to be supplemented by the constraint conditions n^' < 1 and 

The phase diagram of the dipolar system can be easily found by inspecting J^d- 
For iU — E)/t = X oo, the ground state of the system represents a vacuum of 
dipoles. In contrast, for A the ground state is doubly degenerate because 

there are two distinct states with maximal dipole number: {■ ■ ■ d^d^d^ ■ ■ ■)\Q) and 
(• • • d\d\d^^ ■ ■ This immediately suggests the existence of an Ising QCP at some 
intermediate value of A, associated with an order parameter Zi = Y.i{~^Yd'^^di which 
is a density wave of dipoles with a period of 2 lattice spacings. Further analytic 
evidence for an Ising QCP can be obtained by examining the excitation spectra for 
the limiting A regimes, and noting their similarity to those on either side of the critical 
point in the quantum Ising chain yj. 
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Fig. 1.2. Notation as in Fig. II. II (a) A dipole on sites 2 and 3; this state is resonantly coupled 
by an infinitesimal t to the Mott insulator in (a) when E = U. (b) Two dipoles between sites 2 
and 3 and between 4 and 5; this state is connected via multiple resonant transitions to the Mott 
insulator for E = U . 



For A ^ oo, the lowest excited states are sing le dipoles: |^) = t/J|0). There are 
such states (where is the number of sites), and, at A = oo, they are all degenerate 
with energy U — E. The degeneracy is lifted at second order in a perturbation theory 
in 1/X. By a standard approach using canonical transformations, these corrections 
can be described by an effective Hamiltonian, Jfd,eff, which acts entirely within the 
subspace of single dipole states. We find that 



jr,,eff - (f/ - £) E [l^) (^1 + {\i){i\ + \e){i + i\ + \e+i){e\) 



(1.3) 



Notice that, quite remarkably, a local dipole hopping term has appeared in the ef- 
fective Hamiltonian. The constraints (n^ < 1 and n^n^_|_i — 0) played a crucial role 
in the derivation of ( 11.3b . Upon considering perturbations to \£) from the first term 
in jl.H . it initially seems possible to obtain an effective matrix element between 
any two states \£) and |/). However this connection can generally happen via two 
possible intermediate states, \£) djd1,,\0} \i') and \£) |0) — > {£'), and the con- 
tributions of the two processes exactly cancel each other for most £, £' . Only when 
the constraints block the first of these processes is a residual matrix element possi- 
ble. It is a simple matter to diagonalize ^/ eff by going to momentum space; we then 
find a single band of dipole states. The lowest energy dipole state has momentum n: 
the softening of this state upon reducing A is then consistent with the appearance of 
a density wave order with period 2. The higher excited states at large X consist of 
multiparticle continua of this band of dipole states, just as in the Ising chain [IJ. A 
related analysis can be carried out for X — > and the results are similar to those 
for the ordered state in the quantum Ising chain f\\. The lowest excited states form 
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a single band of domain walls between the two filled dipole states, and above them 
are the corresponding multiparticle continua. At an intermediate critical electric field 
Ec = U +1.310t -y/noCno + 1)^ the system undergoes a quantum phase transition lying 
in the Ising universality class 1 17|. 

Having obtained the equilibrium phase diagram for the model, we now consider 
the quench dynamics of the dipoles when the value of the electric field is suddenly 
quenched [4|. We assume that the atoms in the ID lattice are initially in the ground 
state I'J'g) of the dipole Hamiltonian dl.lt with E = Ei ^ Ec- This ground state cor- 
responds to a dipole vacuum. Consider shifting the center of the magnetic trap so that 
the new potential gradient is Ef. If this change is done suddenly, the system initially 
remains in the old ground state. The state of the system at time t is therefore given 
by 

|f(f)) =£c„exp(-/e„f)|«), (1.4) 

n 

where \n) denotes the complete set of energy eigenstates of the Hamiltonian 
with E = Ef, e„ = {n\J^j[Ef]\n) is the energy eigenvalue corresponding to state 
\n), and c„ = (n|f (f = 0)) = {nl^c) denotes the overlap of the old ground state 
with the state |n). (We have set h = 1). Notice that the state |f (f)) is no longer the 
ground state of the new Hamiltonian. Furthermore, in the absence of any dissipative 
mechanism, which is the case for ultracold atoms in optical lattices, |f (f )) will never 
reach the ground state of the new Hamiltonian. Rather, in general, we expect the 
system to thermalize at long enough times, so that the correlations are similar to 
those of Hio[Ef] at some finite temperature. 

We are now in a position to study the dynamics of the Ising density wave order 
parameter 

= ^{^\A\^), (1.5) 
where is the number of sites. The time evolution of O is given by 

Oit) = ^ L cos {m\Y^{~lYdlde\n). (1.6) 

Eq. (11.61 ) is solved numerically using exact diagonalization to obtain the eigenstates 
and eigenvalues of the Hamiltonian Hi£,[Ef]. Before resorting to numerics, it is use- 
ful to discuss the behavior of (9(f) qualitatively. We note that if Ef is close to the 
old ground state will have a large overlap with the new one, i.e., c,„ ~ 5m\ ■ Hence in 
this case we expect (9(f) to have small oscillations about (9(f — 0). On the other hand, 
if Ef ^ Ec, the two ground states will have very little overlap, and we again expect 
(9(f) to have a small oscillation amplitude. This situation is in stark contrast with the 
adiabatic turning on of the potential gradient, where the systems always remain in 
the ground state of the new Hamiltonian Hio[Ef], and therefore has a maximal value 
of ((9) for Ef ^ Ec- In between, for Ef ~ Ec, the ground state | ^) has a finite overlap 
with many states \m), and hence we expect (9(f) to display significant oscillations. 
Furthermore, if the symmetry between the two Ising ordered states is broken slightly 
(as is the case in our studies below), the time averaged value of (9(f) will be non-zero. 
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This qualitative discussion is supported by numerical calculations for finite size 
systems with size = 9, 1 1 , 13. For numerical computations with finite systems, we 
choose systems with an odd number of sites and open boundary conditions, so that 
dipole formation on odd sites is favored, thus breaking the Z2 symmetry. The results 
are shown in Figs. II. 31 - 11. 61 




Fig. 1.3. Evolution of the Ising order parameter in l ll.5t under the Hamiltonian Hio[Ef] for 
no = 1. The initial state is the ground state of HidIEi]. All the plots in this section have U = 40, 
t = I, and Ei = 32, and consequently the equilibrium QCP is at Ec = 41.85. 
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Fig. 1.5. The curve labeled 'dynamic' is the long-time limit {0)t of the Ising order parameter 
in l |1.6t as a function of Er (for N = 11), with other parameters as in Fig. 1 1.1 1 This long-time 



limit can be obtained simply by setting m = n in l ll.6t . For comparison, in the curve labeled 
'adiabatic', we show the expectation value of the Ising order O in the ground state of Hi£)[Ef]; 
such an order would be observed if the value of E was changed adiabatically. Note that the 
dynamic curve has its maximal value near (but not exactly at) the equilibrium QCP Ec = 41 .85, 
where the system is able to respond most easily to the change in the value of E; this dynamic 
curve is our theory of the 'resonant' response in the experiments of Ref. (9) discussed in Sect. 
11.11 In contrast, the adiabatic result increases monotonically with Ef into the E > Ec phase 
where the Ising symmetry is spontaneously broken. 
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Fig. 1.6. Size dependence of the 'dynamic' results in Fig. ll.5l The sizes range from N = 7 to 
N = 15 (as labeled), with the intermediate values A' = 9, 1 1 , 13: (O), decreases monotonically 
with A^. 
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Figure [T3] shows the oscillations of the order parameter (9(f) for different values 
of Ef for = 13. In agreement with our qualitative expectations, the oscillations 
have maximum amplitude when Ef/t w 40 is near the critical value E^/t = 41.85. 
For either Ef <C Ec or Ef 3> Ec, the oscillations have a small amplitude around 
0{t = 0). Furthermore, it is only for Ef k, Ec that the time averaged value of (9(f) is 
appreciable. Figure [L4l shows the system size dependence of the time evolution for 
Ef = U = 40f . We find that the oscillations remain visible as we go to higher sys- 
tem sizes, although they do weaken somewhat. More significantly, the time averaged 
value of (9(f) remains non-zero, and has a weaker decrease with system size. In Fig. 
11.51 we plot the long-time limit of the Ising order parameter, ((9),, as a function of 
Ef, and compare it with (9ad, the value of the order parameter when E reaches Ef 
adiabatically and the wave function is that of the ground state at E = Ef. We find 
that {0)i stays close to (9ad as long as there is a large overlap between the old and 
the new ground states. However, as we approach the adiabatic phase transition point, 
this overlap decreases and {0)t cannot follow (9 ad any more. The deviation of ((9)r is 
therefore a signature that the system is now in a different phase for the new value of 
the electric field. The 'dynamic' curve in Fig. [T3] shows that the Mott insulator has 
a resonantly strong response to an electric field E induced by the proximity of a 
QCP. 

We comment briefly on the nature of the thermodynamic limit, N ^ °° for the 
results in Figs. I1.2l and ll.3l For O^d it is clear that there is a non-zero limit only for 
E > Ec, when it equals the order parameter of the spontaneously broken Ising sym- 
metry. If we assume that the system thermalizes at long times for the dynamic case, 
then {0)t corresponds to the expectation value of the equilibrium order parameter 
in Hij)[Ef] at some finite temperature. In one dimension, it is not possible to break 
a discrete symmetry at finite temperatures, and so the thermodynamic limit of the 
order parameter must always vanish. By this reasoning, we expect {0)t to also van- 
ish in the thermodynamic limit. This is consistent with the results in Fig. 11.61 where 
we show the dependence of the long-time limit ((9),. Our data are at present not 
extensive enough to definitely characterize the dependence of ((9), on A^. 

1.2.2 Infinite range Ising model in a transverse field 

The analysis of the quench dynamics of ID ultracold atoms do not permit an analyti- 
cal description of the long-time value of the order parameter. In particular, the system 
size dependence of the peak height of ((9), is not easy to understand analytically in 
this model. For this purpose, we now consider a simple model system, the infinite 
range ferromagnetic spin- 1/2 Ising model in a transverse field, and study its quench 
dynamics due to a sudden variation of the transverse field. The model Hamiltonian 
is given by ID 

^= - iHs^s^-rJ^sf, (1.7) 

where S" ~ of/2, a = x,y,z, denote the components of the spin-1/2 operator rep- 
resented by the standard Pauli spin matrices cj". Here we assume that / > (fer- 
romagnetic Ising interaction). This Hamiltonian is invariant under the Z2 symmetry 
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S'J S^, S] — > ~S], and S^j — > —5?. (The Z2 symmetry would not be present if there 
was a longitudinal magnetic field coupling to Y.iS^). We take F > without loss 
of generality since we can always resort to the unitary transformation SJ —SJ, 
S-j —5. and 5? ^ S'^, which flips the sign of F but leaves J unchanged. Eq. (11.7b 
can be written as 

H = -^is;,, f - rs-i,, (1.8) 

where Sf,, ^ S^,, = 5f, (1.9) 

and we have dropped a constant {J /2N)Y.i{Sj)^ = J/S from the Hamiltonian in Eq. 
(11.8b . This model has been studied extensively, particularly from the point of view 
of quantum entanglement |18|. Note that this model differs from the one studied in 
Ref. IIT9I . where the spins were taken to be living on two sub-lattices, with Ising 
interactions only between spins on different sub-lattices. 

We begin with a mean field analysis of the thermodynamics of the model de- 
scribed by Eq. ( 11.7b . Denoting the mean field value m = ^, (S'j) /N, the Hamiltonian 
governing any one of the spins is given by 

h = -JmSl, - r5f„, . (1.10) 

This is a two-state problem whose partition function can be found at any temperature 
r. If j3 = 1/ [ksT), we find that m must satisfy the self-consistent equation 



This always has the trivial solution m = 0. In the limit of zero temperature, there is a 
non-trivial solution if F < J /2, with \m\ = (1 /2) ^1 —4r^/J^; the energy gap in that 
case is given by J/2. If F > J/2, we have m = and the gap is given by F — J/2. 
Hence there is a zero temperature phase transition at 77 = 7/2. The Z2 symmetry 
mentioned after Eq. ( 11.7b is spontaneously broken and < 5? > becomes non-zero 
when one crosses from the paramagnetic phase at F > J/2 into the ferromagnetic 
phase r < J/2. 

In the plane of {ksT /J,r /J), there is a ferromagnetic (FM) region in which the 
solution with m 7^ has a lower free energy (the Z2 symmetry is broken), and a para- 
magnetic (PM) region in which m = 0. The boundary between the two is obtained by 
taking the limit m ^ in ( 11.111 ). This gives 2r/J = tanh(j3F/2), i.e.. 



ksT _ r 
J J 



In 



i+2r// \ 
i-2r/yj 



(1.12) 



The mean-field phase diagram is shown in Fig. II. 71 We note that the exact excitation 
spectrum of this model can also be obtained analytically |8 1. 

Having obtained the phase diagram, we now study the quench dynamics across 
the QCP. To begin with, we study the dynamics of the equal-time order parameter 
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Fig. 1.7. Phase diagram of the model in mean field theory. FM and PM denote ferromagnetic 
and paramagnetic regions respectively. 

coiTelation function (EOC) (defined as ((•Sj„,)^) /S^) by changing the transverse field 
r from an initial value I] to a final value Ff suddenly, so that the ground state 
of the system has no time to change during the quench. In this case, just after the 
quench, the ground state of the system can be expressed, in terms of the eigenstates 
\n) of the new Hamiltonian = —{J/^S){S',„,)'^ — FfSf^f as 

IV/) - J^c„|«), (1.13) 

n 

where c„ denotes the overlap of the eigenstate \n) with the old ground state As 
the state of the system evolves, it is given at time t by 

|V/(f)) -E^"^"'"""'!"), (1-14) 

n 

where E„ — {n\,ytff \n) are the energy eigenvalues of the Hamiltonian The EOC 
can thus be written as 

{wmSi,f/S^W{t)) = 1^ c„c„, cos [(£„-£,„)?] {m\{Si,f/S^\n). (1.15) 

m.n 

Eq. ( 11.15b can be solved numerically to obtain the time evolution of the EOC. We 
note that, similar to the case of the dipole model discussed in Sect. 11.2.11 we expect 
the amplitude of oscillations to be maximum when Ff is near f^-. This is verified in 
Fig. 11.81 Here, we have quenched the transverse fields to Ff/J — 0.9,0.01, and 0.4 
starting from Fj/J ~ 2.0. The oscillation amplitudes of the EOC for S = 100, as 
shown in Fig. II. 81 are small for Ff = 0.9 and 0.01, whereas it is substantially larger 
for Ff = 0.4. 
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Fig. 1.8. Dynamics of {{Sfg/)^) /S^ for 5 = 100 after quenching the transverse field to different 
values Ff/J from an initial field i7/y = 2. The oscillation amplitudes are small, as seen from 
the solid (red) and dotted (blue) curves corresponding to iy-// = 0.9 and 0.01 respectively, far 
away from the critical point Fc/J = 0.5. The oscillation is large in the ordered phase near the 
critical point as seen from the dashed (black) curve Ff/J = 0.4. 

Next, to understand the dynamics of the EOC in a little more detail, we study its 
long-time averaged value given by 

0= \im{{{SjJ^{t)))^/S^ 

= ^Lc2(«|(5L)'|«> (1.16) 

for different Ff. Note that the long-time average depends on the product of the over- 
lap of the state |n) with the old ground state and the expectation of in that 
state. From our earlier discussion in Sect. 11.2.11 we therefore expect O to have a 
peak somewhere near the critical point where such an overlap is maximized. This is 
verified by explicit numerical computation of Eq. (|1.16t in Fig. |1.9| for several values 
of S and Fj/J = 2. We find that O peaks around Ff/J = 0.25, and the peak height 
decreases slowly with increasing S. 

To understand the position and the system size dependence of the peak in O, we 
now look at the thermodynamic (large system size) limit; in the present model, this 
is also the large S and therefore classical limit. With this observation, we study the 
classical equations of motion for S — S (cos sin , sin sin , cos 0) for F = Ff. In 
the present model, 5 is a constant. Thus in the classical limit, we need to study the 
equations of motion for 9 and (j). To this end, we note that the classical Lagrangian 
can be written in terms of 6 and as ll20l 
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Fig. 1.9. Plot of the long-time average O as a function of Ff/J for different S. The solid (blue), 
dotted (black), dash-dotted (green) and the dashed (red) lines represent respectively the results 
for 5 = 50, 5 = 100, S = 200 and 5 = 500. O peaks around Fj/J = 0.25, and the peak value 
decreases with increasing S. We have chosen Fi/J = 2 for all the plots. 



L = -5 [1 -COS0] ^ - ^[0,0]. (1.17) 
This gives the equations of motion 

r ■ A 

^=^^^^"'^' 

d<i) J 

— = cos + Ff cote cos . (1.18) 

at 2 

Eq. (11.181) has to be supplemented with the initial condition that S^i,t ~ S ai t =Q. 
The condition 5f^_„ S corresponds to0 = ;r/2,^=O which is also a fixed point 
of (11.18b . Therefore we shall start from an initial condition which is very close to 
the fixed point: = n:/2 — £, ^ = e, where e is an arbitrarily small constant. Further, 
since the motion occurs on a constant energy surface after the quench has taken place, 
we have 

F/- = - cos^e + Ff sin0cos0 . (1.19) 
Using (II.I8I 1 and (|1.19l l. we get an equation of motion for Q in closed form. 



dd 



Ffs.i\\^{e)- [r^-|cos2 



at sin 

It can be seen that the motion of Q is oscillatory and has classical turning points at 
0, =sin"' (|l -4r/-/7|) and 02 = 7r/2. One can now obtain {{Sl„f)^ = (cos2 0)^ 
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from ( foot . 



(cos^e)^ 



(1.21) 



where .yV 




cos^ 



(1.22) 



1 



and ^ 



(1.23) 



When trying to evaluate we find that the integral has an end-point singularity at 
02; this can be regulated by a cut-off 77 so that 62 = — fj. With this regularization, 



^ — —J\s\(r\)l y Tf (j — IFf^ /2. The cut-off used here has a physical meaning and 

is not arbitrary. To see this, note that the angles (0,0) define the surface of a unit 
sphere of area 4-7t. This surface, for a system with spin S, is also the phase space 
which has 25+1 quantum mechanical states. For large S, the area of the surface 
occupied by each quantum mechanical state is therefore 4;r/(25'+ 1) ~ In/S. In 
other words, each quantum mechanical state will have a linear dimension of order 
1 / V5; this is how close we can get to a given point on the surface of the sphere. 
Note that this closeness is determined purely by quantum fluctuation and vanishes 
for 5 — > 00. Thus rj, which is also a measure of how close to the point 9 = k/2 we 
can get, must be of the order of 1 / \/S; this determines the system-size dependence 
of (cos^ 0)^. Using (II.2II 1. we finally get 



Eq. (|1.24| i is one of the main results of this section. It demonstrates that the long- 
time average of the EOC must be peaked at Ff/J = 0.25 which agrees perfectly with 
the exact quantum mechanical numerical analysis leading to Fig. 11.91 Moreover, it 
provides an analytical understanding of the S (and hence system size) dependence of 
the peak values of we conclude that the peak in O vanishes logarithmically with 
the system size S. Such a slow variation with S shows that it might be experimentally 
possible to observe an experimental signature of a QCP for a possible realization of 
this model with ultracold atoms where ^ 1 0^ — 10^ 1 8 1 . 

The results obtained in this section can be tested in two kinds of experimental 
systems. One class of systems are those with long-range dipole-dipole interactions 
such as KH2PO4 or Dy(C2H5S04)39H20 (2H which exhibit order-disorder transi- 
tions driven by tunneling fields. The other class of systems are two-component Bose- 
Einstein condensates where the inter-species interaction is strong compared to the 
intra-species interaction; the relative strengths of these interactions can be changed 
by tuning the system to be near a Feshbach resonance as discussed for the ^^K—^^Rb 
system in Refs. Il22ll23]l24ll25l . The quench dynamics that we have discussed can 
be realized by applying a radio frequency pulse to the system and suddenly changing 
the frequency of the pulse. 

We end this section with the observation that the resonant response of the order 
parameter during quench dynamics has been found in two very disparate models 




(cos' 





(1.24) 



jnn{S) 
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(ID ultracold atoms and infinite dimensional Ising ferromagnet) and therefore seems 
to be a universal signature of the QCP through which the system passes during its 
evolution. 



1.3 Non-adiabatic dynamics 

In recent years, there have been extensive studies of what happens when a pa- 
rameter A in the Hamiltonian of a quantum system is varied in time slowly (non- 
adiabatically) so as to take the system through a QCP. A quantum phase transition is 
necessarily accompanied by diverging length and and time scales, or, equivalently, a 
vanishing energy gap between the ground state and the first excited state [ 1 1. A con- 
sequence of this is that the system fails to be in the adiabatic limit when it crosses 
a critical point. Namely, when X is varied across the QCP located at A = at a 
finite rate given by 1 /t (where T will be called the quench time), the system fails to 
follow the instantaneous ground state in a finite region around A^. As a result, de- 
fects are produced ifTOl l26l . For a slow quench (for instance, for values of T much 
larger than the inverse band width) which takes the system across a QCP in a linear 
way, it is well-known that the density of defects n scales as a power of the quench 
time, n ~ where v and z are respectively the correlation length and 

the dynamical critical exponents characterizing the critical point [12^, "271. A theo- 
retical study of a quench dynamics requires a knowledge of the excited states of the 
system. Hence, early studies of the quench problem were mostly restricted to quan- 
tum phase transitions in exactly solvable models such as the ID Ising model in a 
transverse field [4, 28 29], the ID XY spin-1/2 model ll30l[33]| . quantum spin chains 
EHmEl, the Bose-Hubbard model |E21, the Fahcov-Kimball model iMl and ID 
spinless fermionic chains f37l. Experimentally, trapped ultracold atoms in optical 
lattices provide possibilities of realization of many of the above-mentioned systems 
Il38ll39l . Experimental studies of defect production due to quenching of the magnetic 
field in a spin-1 Bose condensate have also been undertaken [40 1. 

A class of models in which the above power law scaling can be derived easily 
is one in which, due to the existence of a mapping to a system of non-interacting 
fermions, the system decomposes into a product of two-level systems. For instance, 
this occurs in the ID XY spin-1/2 model and in the two-dimensional (2D) spin-1/2 
Kitaev model. In both these cases, it turns out that the Hamiltonian is given by a sum 
of terms of the form 

//k = a(k)(4ck + cl^ck) + 4*(k)4ct^ + 4(k)c_kCk, (1.25) 

where k runs over halfihe Brillouin zone (BZ), and a(k) is real. This Hamiltonian 
acts on a space spanned by four states, namely, the empty state |0), two one-fermion 
states c^|0) = |k) and c^i^jO) = | — k), and a two-fermion state c^c^i^jO) = |k, — k). 
Both the one-particle states are eigenstates of //k with the same eigenvalue a(k). On 
the other hand, the states |0) and |k, — k) are governed by a 2 x 2 Hamiltonian given 
by 
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'-(..;)2i,'S)- 

The eigenvalues of this are given by o:(k) ± ^ a?-{\&) + |zi(k)p; since the lower 
eigenvalue is less than a(k), the ground state lies within the subspace of the states 
|0) and |k, — k) . We observe that the Hamiltonian in ( 11.25b does not mix the states |0) 
and |k, — k) with the states |k) and | — k), even !7a(k) and 4(k) change with time. 
Hence, if we start at time f ^ — oo with a linear superposition of |0) and |k,— k), 
we will end at f oo with a superposition of the same two states. In that case, it is 
sufficient to restrict our attention to the Hamiltonian for a two-level system given in 
(11.261 ). We can rewrite that as 

K = a(k) h - a(k) al + 4(k) (7+ + 4*(k) , (1.27) 

where is the identity matrix, and C7^ and C7* = (a' ±i(7^)/2 denote the Pauli 
matrices. We can ignore the term a(k)/2 in ( ll.27l i since this only affects the wave 
function by a time-dependent phase factor. 

Let us now consider what happens if a(k) varies linearly with time. Then the 
total Hamiltonian //j for all the two-level systems can be written as 

Hd = ^ hk, where hk = - e(k) al + A{k) C7+ + A*{k) a^, (1.28) 

k ^ 

where d is the number of spatial dimensions, and the sum over k runs over half 
the BZ. If e(k) > 0, the ground state is given by |0) as f ^ — oo and by |k, — k) as 
? ^ oo. If we begin with the state |0) at f = —°a and evolve the system using the time- 
dependent Schrodinger equation, we end at f = in a state which is a superposition 
of states |0) and |k,— k) with probabilities pk and 1 — pk, where pk is given by the 
Landau-Zener expression ||4T1 

Pk = e-'^^I^WI'/^W. (1.29) 

Note that pk 1 for t ^ (sudden quench) and for t ^ °o (adiabatic quench), 
as expected. Let us now assume that the gap function A (k) vanishes at some point 
ko in the BZ as |4(k)| ~ flo|k — ko|', while e(ko) — bo is finite; this corresponds 
to a QCP with an arbitrary value of z, but with zv = 1. The density of defects n in 
the final state is given by the density of fermions n — J d'^kpk- In the adiabatic Umit 
T ^ oo, this is given by 

/ e-''™olk-koP7''o _ i/t^/(2z)^ (1.30) 

Jk~ko 

which is the expected result for zv = I. 

It is useful to note that the derivation of the scaling law in ( 11.301 ) does not require 
a knowledge of the precise functional form given in ( 11.29b . It is enough to know that 
Pk must be a function of the form f{T\A (k)p/£(k)), where f{x) 1 for T ^ and 
— + for T oo; these limiting values follow from general properties of the time- 
dependent Schrodinger equation. The argument of the function / can be derived by 
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considering the equation id\i/k{t)/dt — /!kVk(0' performing a phase re-definition 
to change A{k) to |4(k)|, multiplying both sides by ^ T/£(k) and rescahng t to 
f-\/T/e(k). This effectively converts /zu to the form ta^ + |zi(k)|-\/T/e(k)(c7|^ + 

CT]^); hence the probability pk of starting in the ground state atf = — oo and 

ending in the same state (which is the excited state) at t = °o must be a function of 
T|4(k)|Ve(k). 

We can generalize the above results to a QCP with arbitrary values of z and v. 
We consider a generic time-dependent Hamiltonian H{t) = H[X{t)], whose states 
are labeled by |k), and |0) denotes the ground state. If there is a second order phase 
transition, the basis states change continuously with time during this evolution and 
can be written as \yf{t)) = Lk'2k(0l'^['^(0])- The defect density can be obtained in 
terms of the coefficients flk(0 as n = Lk/o Wk{t — * hence one gets l,12j 



d''k 



dX (k|-^|0) 

3 UA 



(1.31) 



where 5a)k(A) = a)k(^) — ftJD(A) are the instantaneous excitation energies. Follow- 
ing Ref. lfT2l . we note that near a QCP, 5a)k(A) = i4F(/i/|k| '), where A is the energy 
gap, z is the dynamical critical exponent and F(x) ~ 1 /x for large x\ we have assumed 
here that the gap vanishes as k = 0. Also, since the quench term vanishes at the criti- 
cal point as A - |A|''^, one can write 5c(Jk(A) ^ lAh'^FdAh'^^/lkl^') where - l/x 
for large x. Further, one has (k|^|0) — |k|^^'G(4/|k|^') near a critical point where 
G(0) is a constant. This allows us to write (k| ^ |0) = A^'^"^ |k|-^G'(A^^/|k|') where 
G'(0) is a constant ([T][l2l. Substituting these in (11.311 1 and changing the integration 
variables to Tj = W(^^+i)|k| and t, = |k|-i/v;t, we find that 

We will now discuss two major extensions of the above results: (i) what happens 
if the system is taken across ad — m dimensional quantum critical surface instead of a 
QCP |fT3l , and (ii) what happens if the quenching across a QCP is non-linear in time 
lfT4ll42l . We will show that in both cases, the defect density still scales as a power 
of the quench time, but the power is not equal to the universal value dv / {zv + 1 ) 
mentioned above; rather it depends on other parameters such as m or the degree of 
non-linearity. 



1.3.1 Quenching across a critical surface 

When a quench takes a quantum system across a critical surface rather than a critical 
point, the density of defects scales in a different way with the quench time. To give a 
simple argument, consider a cZ-dimensional model with z = v = 1 which is described 
by the Hamiltonian given in ( 11.281 ). Suppose that a quench takes the system through a 
critical surface ofd — m dimensions. The defect density for a sufficiently slow quench 
is then given by ||4T| n ~ J^^^^^^ '^^'^^^^' where /(k) — |zi(k)p/|£(k)| vanishes on 
the d — m dimensional critical surface. We can then write 
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d'kexpl -TTT £ gapkakp] - 1/W, (1.33) 
a.p=l 

where a,j3 denote one of the m directions orthogonal to the critical surface, and 
8aj3 = [5 V(k)/<5^a<3^j3]kGcriticai surface- Note that this result depends only on the 
property that /(k) vanishes on a d — m dimensional surface, and not on the precise 
form of /(k). For general values of z and v, we note that the Landau-Zener type of 
scaUng argument yields A ^ l / t''*'/(^*'+ ' ' , where A is the energy gap 1 1 2 1 . When one 
crosses a d — m dimensional critical surface during the quench, the available phase 
space Q for defect production scales as i2 k'" ^ A'"l~ ^ this leads 

to n ^ 1/ -j-'n*'/(2*'+i) For a quench through a critical point where m — d, we retrieve 
the results ofRef. |T2l. 

To give an example of a quench across a critical line, let us consider a model 
which was proposed recently by Kitaev. This is a 2D spin- 1/2 model on a honeycomb 
lattice as shown in Fig. ll.lOl the Hamiltonian is given by B6l 

Hk= £ + -^zfTj^i./CTj, + ho]ja]j^, ), (1.34) 

j+/=even 

where j and / denote the column and row indices of the honeycomb lattice. This 
model has been studied extensively and it exhibits several interesting features ll47l 
l48l |49l l50l |5T] |52]. It provides a rare example of a 2D model which can be exactly 
solved using a Jordan- Wigner transformation |46]|47]|5Tl|52l. It has been shown in 
Ref. P6l that the presence of magnetic field, which induces a gap in the 2D gapless 
phase, leads to non-Abelian statistics of the low-lying excitations of the model; these 
excitations can be viewed as robust qubits in a quantum computer 1531 . 

The Jordan- Wigner transformation of the Kitaev model to a model of non- 
interacting fermions works as follows. One can write 



fly/ = 



n svenj + l, 



^j' = ri ^J/ for odd; + /. (1.35) 

where the Uji and bji are Majorana fermion operators (and hence Hermitian) obeying 
the anticommutation relations {fly/,fly//} = {bji,bjiii} = djfdii'. This transformation 
maps the spin Hamiltonian in (ll.34l i to a fermionic Hamiltonian given by 

^nfln-Mi + h ^nfln+M2 + J^Dn bnOn], (1.36) 

n 

where n = %/3/ ni -t- (-^/ + jj) «2 denote the midpoints of the vertical bonds. Here 
ni,n2 run over all integers so that the vectors n form a triangular lattice whose ver- 
tices lie at the centers of the vertical bonds of the underlying honeycomb lattice; the 
Majorana fermions a„ and ba sit at the top and bottom sites respectively of the bond 
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Fig. 1.10. Schematic representation of the Kitaev model on a honeycomb lattice showing the 
bonds 71, J2 and Jt,. Schematic pictures of the ground states, which correspond to pairs of 
spins on vertical bonds locked parallel (antiparallel) to each other in the limit of large negative 
(positive) /3, are shown at one bond on the left (right) edge respectively. Mj and M2 are 
spaiming vectors of the lattice, and a and b represent inequivalent sites. 



labeled n. The vectors Mi = and M2 — jj are spanning vectors for 

the reciprocal lattice, and Da can take the values ± 1 independently for each n. The 
crucial point that makes the solution of Kitaev model feasible is that Dn commutes 
with Hk, so that all the eigenstates of Hk can be labeled by a specific set of values 
of Dn. It has been shown that for any value of the parameters 7,, the ground state 
of the model always corresponds to D„ = 1 on all the bonds. Since D„ is a constant 
of motion, the dynamics of the model starting from any ground state never takes the 
system outside the manifold of states with Dn = 1 . 

For Dn = 1, it is straightforward to diagonalize Hk in momentum space. We 
define Fourier transforms of the Majorana operators an as 




(1.37) 



where is the number of sites (assumed to be even, so that the number of unit 
cells N/2 is an integer), and the sum over k extends over half the Brillouin zone 
of the honeycomb lattice. We have the anticommutation relations {flk,fl^/} = Sf^k', 
{ak,a^i} — 0, and similarly for b^, and b^. We then obtain Hk = Y,k Wk^kVk, where 
Vk ~ (''k'^k)' ^k can be expressed in terms of Pauli matrices a^-^-^ as 
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/zk = 2 [Ji sin(k • Ml ) - J2 sin(k • M2)] 

+ 2 [73 + /icos(k-Mi) + 72Cos(k-M2)] a^. (1.38) 

The energy spectrum of Hk consists of two bands with energies 

^±2 [(7isin(k-Mi)-y2sin(k-M2)f 

+ (/3 +^1 cos(k • Ml ) +y2Cos(k • M2))2] '/^ (1.39) 

We note that for |7i —J2\ <J3<Ji +J2, these bands touch each other so that the 
energy gap — — E^ vanishes for special values of k leading to a gapless phase 
of the model ESI El HI] ED ■ 

We will now quench 73(f) —Jt/T at a fixed rate 1/t, from —00 to 00, keeping J, 
J I and J2 fixed at some non-zero values; we have introduced the quantity J to fix the 
scale of energy. We note that the ground states of Hk corresponding to 73 — > —00(00) 
are gapped and have ctj^ctj^^j = 1(— 1) for all lattice sites (y,/). To study the state 
of the system after the quench, we first note that after an unitary transformation 
U — exp(— /(7i;r/4), one can write Hk = Lk W^^'kVk' where h'-^ ~ Uh^U^ is given 
by 

/z'k = 2 [7isin(k-Mi) - 72sin(k-M2)] 

+ 2 [y3(f)+yicos(k-Mi)+72Cos(k-M2)] (7^ 

(1.40) 

Hence the off-diagonal elements of h'^ remain time-independent, and the problem of 
quench dynamics reduces to a Landau-Zener problem for each k. The defect density 
can then be computed following a standard prescription BTl 



1 



A Jk 

= e-2jiT [yisin(k-Mi)-72sin(k-M2)]V^^ (141) 

where A — 4a:^/(3\/3) denotes the area of half the Brillouin zone over which the 
integration is carried out. Since the integrand in ( 11.411 ) is an even function of k, 
one can extend the region of integration over the full Brillouin zone. This region 
can be chosen to be a rhombus with vertices lying at {kx,ky) = (±27r/\/3,0) and 
(0, ±271/3). Introducing two independent integration variable vi,V2, each with a 
range < vi , V2 < 1, one finds that 

k,^2n^--^ , ky^2n^^. (1.42) 

Such a substitution covers the rhombus uniformly and facilitates the numerical inte- 
gration necessary for computing n. 

A plot of n as a function of the quench time Jx and a = tan^' (72/^1) (we have 
taken /jp] = 7cos a[sin a]) is shown in Fig. ll.l ll We note that the density of defects 
produced is maximum when J\ = J2. This is due to the fact that the length of the 



1 Non-equilibrium dynamics of quantum systems 



21 



gapless line through which the system passes during the quench is maximum at this 
point. This allows the system to remain in the non-adiabatic state for the maximum 
time during the quench, leading to the maximum density of defects. For J\/Jt, > 
2J2/J3, the system does not pass through a gapless phase during the quench, and the 
defect production is exponentially suppressed. 




The density of defects is maximum at Ji = J2- 

For sufficiently slow quench InJT ^ 1, /7k is exponentially small for all values 
of k except in the region near the line 

Ji sin(k-Mi) - 72sin(k-M2) = 0, (1.43) 

and the contribution to the momentum integral in (11.411 1 comes from values of k 
close to this line of zeroes. We note that the line of zeroes where = 1 precisely 
corresponds to the zeroes of the energy gap Zik as ^3 is varied for a fixed J2 and 
Ji . Thus the system becomes non-adiabatic when it passes through the intermediate 
gapless phase in the interval |7i —J2\< Ji{t) < J\ +-/2- It is then easy to see, by 
expanding about this line that in the limit of slow quench, the defect density scales 
as n '-^ 1 / We thus see that the scaling of the defect density with the quench rate 
when the system passes through a critical line in momentum space is different from 
the situation where the quench takes the system through a critical point. The Kitaev 
model is an example of a system in which d = 2, m = I, and z— V ~ 1; this gives 
rise to a defect density scaling as 1 /Vt- 

Before ending this section, it is interesting to consider another aspect of the sys- 
tem after the quench. Since the time evolution of the system is unitary, it will al- 
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ways be in a pure state. However, for each value of k, the wave function is given by 
VI -/7kV^2ke"'^2'' + y/PkWike^'^^''' , where En, (£2^) = °° (-°°)- As a result, the 
final density matrix of the system will have off-diagonal terms involving v4k V^ik ^"'^ 
VikVik which vary extremely rapidly with time; their effects on physical quantities 
will therefore average to zero. Hence, for each momentum k, the final density matrix 
Pk is effectively diagonal like that of a mixed state t30il . where the diagonal entries 
are time-independent ast^°o and are given by \ — pk and p^. Such a density matrix 
is associated with an entropy which we will now calculate. The density matrix of the 
entire system takes the product form p = (^Pk- The von Neumann entropy density 
corresponding to this state is given by 

.V = - 1 j d\[{\-pM^-P\d + /'kln/^k], (1.44) 

where the integral again goes half the Brillouin zone. Let us now consider the depen- 
dence of this quantity on the quench time T f3T|. If T is very small, the system stays 
in its initial state and pk will be close to 1 for all values of k; for the same reason, 
(Oq) will remain close to 1 . If t is very large, the system makes a transition to the 
final ground state for all momentum except near the line described in (11.43b . Hence 
Pk will be close to for all k except near that line, and ((9o) will be close to -1. In 
both these cases, the entropy density will be small. We therefore expect that there 
will be an intermediate region of values of T in which s will show a maximum and 
(Oo) will show a crossover from — 1 to 1 . A plot of s as a function of Jz and a shown 
in Fig. I I.I2I confirms this expectation. We find that the entropy reaches a maximum 
for an intermediate value of Jt where (Co) crosses over from — 1 to 1 for all values 
of a. 

1.3.2 Non-linear quenching across a critical point 

Let us now consider what happens if we start with a Hamiltonian similar to the one 
given in (11.281 1, except that 

/!k(f) = {X{t) + b{\.))al + 4(k) (j+ + 4*(k) , (1.45) 

where A (f) = Ao|f/T|"sign(f) is the quench parameter; a~l corresponds to a Unear 
quench. The instantaneous energies of the Hamiltonian in ( ll.45l l are given by 

£(k) = ± ^(A(f)+Kk))2 + |4(k)|2. (1.46) 

These energy levels touch each other at f = fo and k = ko, so that \A (k)| ^ |k — ko| 
and |fo| = T|Z7(ko)/Ao|^''" = Tg'^", where^= |/7(ko)/Ao| is a non-universal model- 
dependent parameter. At this point the energy levels cross and we have a QCP with 
Z = V = 1. Note that the critical point is reached at f = only if /^(ko) vanishes. 

Let us first consider the case where fe(ko) = so that the system passes through 
the critical point at f = 0. In what follows, we shall assume that |4 (k)| ~ |k — ko| and 
b{\i.) ~ |k — kol'i at the critical point, where zi > 1 so that E ^ |k — ko| and z = 1. 
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0.06 




1.00 



Fig. 1.12. Plot of the entropy density i as a function of Jx and a = tan ' (Ji/Ji)- The entropy 
density peaks when (Oo) crosses from —1 to 1. 



In the rest of the analysis, we will scale f fAo, T tXq, A{k) A{k)/Xo, and 

To obtain the probability of ending in the excited state at f = oo, we study the 
time evolution of the system governed by the Schrodinger equation idxffk/dt = /zkVk- 
This leads to the equations 



icik = (|f/T|"sign(f) + /7(k))cik + A{k) C2k, 
ic2k = -(|f/T|"sign(f)+/7(k)) C2k + A*{k) cik, 

where cik(2k) = '5fCik(2k)- To solve these equations, we define 

JJ'dt'i\t'/T\'',ignit')+b(k)) 



(1.47) 



Cik = cik e 



,-i/'./f'(|r'/T|«sign(r')+/'(k)) 



(1.48) 



C2k = C2k e 

Substituting ( 11.481 ) in ( 1 1 .47l l and eliminating c'ji^ from the resulting equations, we get 
c'lk - 2i [|f/Trsign(f) + fe(k)] c'lk + |4(k)|2 c\k = 0. (1.49) 
Now we scale t so that ( |1.49t becomes 

c'lk - 2i [|r|«sign(f) + /?(k)T«/(«+')] c[k + \A{k)\h^"/'^"+^^ c[k = 0. (1.50) 



From d 1.50b we note that since cik and c[f^ differ only by a phase factor, pk must be 
of the form 
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= /[/'(k)T^,|4(k)|2T^], (1.51) 

where / is a function whose analytical form is not known for a ^ 1. Nevertheless, 
we note that for a slow quench (large t), pk becomes appreciable only when the 
instantaneous energy gap, as obtained from ( 11.46b . becomes small at some point of 
time during the quench. Consequently, / must vanish when either of its arguments 
are large: f{°o,a) — f{a,°o) — for any value of a. Thus for a slow quench (large t), 
the defect density n is given by 

t/'';t/[fe(k)T^,|zi(k)|2T^], (1.52) 

BZ 

and receives its main contribution from values of / near k = ko where both b{k) and 
A (k) vanish. Thus one obtains, after extending the range of the integration to oo, 

(1.53) 



d"kf |k-ko|•■lTS+T;|k-ko|^TS^ 

Now scaling k (k — ko) we find that 

^i''/t/(|k|^iT«(i-''')/(«+i);|k|2) 



/ d^'k f{0;\lif-) - T^s+T, (1.54) 



where, in arriving at the last line, we have used zi > 1 and T — > oo. (If Zi = 1, the 
integral in the first line is independent of T, so the scaling argument still holds). Note 
that for a = 1, Eq. ( I1.54l i reduces to its counterpart for a linear quench [IZj. It turns 
out that the case zi < 1 deserves a detailed discussion which is given in Ref. lfT4l . 

Next we generalize our results to a critical point with arbitrary values of z and v. 
We use arguments similar to those given in the discussion around Eq. ( 11.311 ). namely. 



d''k 



dX (k|-^|0) 



(1.55) 



In the present case, the quench term vanishes at the critical point as 4 ^ |A |"-'^ for a 
nonlinear quench, and we can write 

5cOk(A) - |Ar'-^F(|Ar^7|k-kon, (1.56) 

where /"(jc) - 1 /x for large x Further, (k|^|0) = |k - ko|"-'G(4/|k - koh') near a 
critical point, where G(0) is a constant. This allows us to write 

where G'(0) is a constant lHJIIl]. Substituting ( 11.56b and ( 11.57b in ( 11.551 ) and chang- 
ing the integration variables to 77 = T«^/(«^^+i)|k- ko| and E, = |k - ko|-'/(«^)A, 
we find that 
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Next we consider the case where the quench term does not vanish at the QCP 
for k = ko. We again consider the Hamiltonian h^{t) in ( 11.451 ). but now assume that 
the critical point is reached sit — Iq^Q. This renders our previous scaling argument 
invalid since 4(ko) — but /^(ko) ^ 0. In this situation, |fo/T| — g^l" so that the 
energy gap AE may vanish at the critical point for k — ko. We now note that the 
most important contribution to the defect production comes from times near fo and 
from wave numbers near kQ. Hence we expand the diagonal terms in /!k(r) about 
t — tQ and k = ko to obtain 



(1.59) 



where b'{5k) represents all the terms in the expansion of bik) about k = ko, and we 
have neglected all terms 

R„ = (a-n + l)(a-n + 2)...(a)g(«-"V«|(f-fo)/T|"sign(f)/n! (1.60) 

for « > 1 in the expansion of A (f ) about fo. We shall justify neglecting these higher 
order terms shortly. 

Eq. ( 11.591 ) describes a linear quench of the system with Teff(a) = t/ (ag^"^ ''/"). 
Hence one can use the well-known results of Landau-Zener dynamics BTl to write 
an expression for the defect density, 

«~ / d^'kpk^ [ d''k exp[~7t\A(k)\^T^s{a)]. (1.61) 
Jbz Jbz 

For a slow quench, the contribution to n comes from k near ko; hence 

n - Teff(a)-''/2 ^ fa^(«-i)/«/Ty'/'. (1.62) 



Note that for the special case a = 1, we recover the familiar result n '--^ i^'^l^, and the 
dependence of n on the non-universal constant g vanishes. Also, since the quench is 
effectively linear, we can use the results of Ref. lfT2l to find the scaling of the defect 
density when the critical point at f = fo is characterized by arbitrary z and v. 



(1.63) 



Next we justify neglecting the higher order terms R„. We note that significant 
contributions to n come at times f when the instantaneous energy levels of H in 
( |1.59t for a given k are close to each other, i.e., (f — fo)/T Zi (k). Also, for a slow 
quench, the contribution to the defect density is substantial only when is signifi- 
cant, namely, when \ A (k)p ^1/ Teff(a). Using these arguments, we see that 

Rn/Rn-x = (a-n + l)g-'/"(f-fo)/(nT) - (a-n + l)/(«Vi). (1.64) 
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Thus we find that all higher order terms i , which were neglected in arriving at 
( 11.621 ). are unimportant in the limit of slow quench (large t). 

The scaling relations for the defect density n given by (11.581 1 and ( 11.63b represent 
the central results of this section. For such power law quenches, unlike their linear 
counterpart, n depends crucially on whether or not the quench term vanishes at the 
critical point. For quenches which do not vanish at the critical point, n scales with 
the same exponent as that of a linear quench, but is characterized by a modified non- 
universal effective rate Teff(o:). If, however, the quench term vanishes at the critical 
point, we find that n scales with a novel a-dependent exponent adv/{azv + 1). For 
a = 1, Teff(a) = T and adv/{azv + l) = dv/{zv + l); hence both ( 11.581 ) and ( 11.63b 
reproduce the well-known defect production law for linear quenches as a special case 
lfT2ll . We note that the scaling of n will show a cross-over between the expressions 
given in (11.581 ) and ( 11.631 ) near some value of T = To which can be found by equating 
these two expressions; this yields To ^ |/7(ko)|^'^^'/". For a > 1, the scaling law 
will thus be given by Eq. ( 11.58b (( 11.631 )) for T ^ (>-)to. We also note here that 
the results of this section assumes that the system passes from one gapped phase to 
another through a critical point and do not apply to quenches which take a system 
along a critical Hne lfT3]|43ll44l|45i . 

To illustrate the form of defect scaling for a non-linear quench, let us consider 
the ID spin- 1/2 Kitaev model which is governed by the Hamiltonian 

H = + 725jCi)' (1-65) 

/eeven 

where S" — crf /2. Using the standard Jordan- Wigner transformation, this can be 
mapped to a Hamiltonian of non-interacting fermions 

k 

where h/^ = —2 (y_sinfe T3 +J+cosk T2). (1.66) 

Here J± = 7i ±72, and Xj/k = (^i (fc),C2(fe)) are the fermionic fields. We now per- 
form a quench by keeping fixed and varying the parameter J with time as 
J-{t) = 7|f/T|"sign(f). We then pass through a QCP at f = at the wave num- 
ber k = 7t/2. From Eq. (11.58b we expect the defect density to go as n '--^ ^-a/(a+i) 
since v = z = 1 for this system. To check this prediction, we numerically solve the 
Schrodinger equation id\i/{k,t)/ dt = hii{t)\i/{k,t) and compute the defect density 
n — Jq [dk/n) pk as a function of the quench rate T for different a, with fixed 
J+/J= 1. A plot of ln(n) vs ln(T) for different values of a is shown in Fig. 11.131 
The slopes of these lines, as can be seen from Fig. 11.131 changes from —0.67 towards 
— 1 as a increases from 2 towards larger values. This behavior is consistent with the 
prediction of ( 11.58b . The slopes of these lines show excellent agreement with (11.58b 
as shown in the inset of Fig. 11.131 

To illustrate what happens if the QCP is crossed at a time t which is different 
from 0, we consider the ID Ising model in a transverse magnetic field described by 

^^sing = - -/ (E S]S\^i + g E ^?)' (1-67) 
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Fig. 1.13. (Color online) Plot of ln(n) vs ln(T) for the ID Kitaev model for a = 2 (black solid 
line), a = 4 (red dotted line), a = 6 (blue dashed line) and a = 8 (green dash-dotted line). 
The slopes of these lines agree reasonably with the predicted theoretical values —a/{a + 1) 
as shown in the table. 



where J is the strength of the nearest neighbor interaction, and g ^ li/ J is the di- 
mensionless transverse field. In what follows, we shall quench the transverse field as 
g{t) — |f/T|"sign(f ) and compute the density of the resultant defects. 

We begin by mapping //king to a system of free fermions via the Jordan- Wigner 
transformation 



k 



{g — cosk) al + sin^cj/]. 



(1.68) 



If the field g is varied with time as g{t) ~ gQ\t / T|"sign(f ), the system will go through 
two QCPs at g = 1 and — 1 . The energy gap vanishes at these QCPs at = fco = and 
n. As a result, defects are produced in non-adiabatic regions near these points. For 
this model, the QCP is at f = fo ^ and z = V — I- Hence, Tgff = t/o: for both the 
QCPs. From (11.63b . therefore, we expect the defect density produced in this system 
to be given by n ^ (r/a)"'/^. 

To verify this, we numerically solve the Schrodinger equation id\i/ic{t)/ dt = 
hk{t)Wk{t) and obtain the probability p/^ for the system to be in the excited state. 
Finally, integrating over all k within the Brillouin zone, we obtain the defect density 
n for different values of a > 1 with fixed T. The plot of « as a function of a for 
T = 10, 15 and 20 is shown in Fig. 11.141 A fit to these curves gives the values of the 
exponents to be 0.506 ± 0.006, 0.504 ± 0.004 and 0.505 ± 0.002 for T = 10, 15 and 
20 respectively which are remarkably close to the theoretical value 1 /2. The system- 
atic positive deviation of the exponents from the theoretical value 1 /2 comes from 
the contribution of the higher order terms neglected in the derivation of ( 11.62b and 
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Fig. 1.14. (Color online) Variation of the defect density n with the quench exponent a for 
representative values of T = 10 (black solid line), T = 15 (red dashed line) and T = 20 (blue 
dotted line). A polynomial fit of the form n = aa!' yields exponents which are very close to 
the theoretical result 1/2 for all values of t. 

(11.631 1. We note that the region of validity of our linear expansion, as can be seen 
from Fig. 03 

grows with T which is in accordance with the result in (11.641) . 
1.3.3 Experimental realizations 

The validity of our results can be checked in a variety of experimental systems. We 
first observe that all our results have been obtained at zero temperature with the 
assumption that the system does not relax significantly during the quenching process 
and until the defect density has been measured. This might seem rather restrictive. 
We note however that systems of ultracold atoms in optical or magnetic traps and/or 
optical lattices can easily satisfy the required criteria since they have a very long 
relaxation time which often gets close to the system lifetime ll38l . We will briefly 
list some possible experiments here. First, there has been a proposal for realizing 
the Kitaev model using an optical lattice||39]|. In such a realization, all the couplings 
can be independently tuned using separate microwave radiations. In the proposed 
experiment, one needs to keep 73 = and vary 71(2) — -^(1 =t |r/T|"sign(f))/2, so 
that y+ remains constant while J varies in time. The variation of the defect density, 
which in the experimental set-up would correspond to the bosons being in the wrong 
spin state, would then show the theoretically predicted power law behavior in ( 11.581 ). 
Secondly, a similar quench experiment can be carried out with spin-1 bosons in a 
magnetic field described by an effective Hamiltonian //gff = C2no(S)^ +ciB^(5'^) 
P-OI . where C2 < and no is the boson density. Such a system undergoes a quantum 
phase transition from a ferromagnetic state to a polar condensate at B* = \/|c2|no/ci. 
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A quench of the magnetic field — BqIi/t]" would thus lead to a scaHng of the 
defect density with an effective rate Teff(a) = T/(ag'"^^''''"), where g — |c2|no/ci. 
A measurement of the dependence of the defect density « on a should therefore 
serve as a test of the prediction in ( 11.631 1. Finally, spin gap dimer compounds such 
as BaCuSi206 are known to undergo a singlet-triplet quantum phase transition of 
the Bose-Einstein condensation type at Be ~ 23T; the critical exponents for this are 
given by z = 2 and V — 2/d. Experimentally, the exponent v appears to be 2/3 above 
a temperature window of 0.65K to 0.9K ||54l and 1 below that temperature window 
due to a dimensional reduction from d = 3to d = 2 |[55l . Thus a nonlinear quench of 
the magnetic field through its critical value B = Be +Bo|f/T|"sign(f) should lead to 
a scaling of the defects n - x-(><^l{'^»+^) for ^ 3, V 2/3, and « - T-2a/(2a+i) foj. 

= 2, V = 1 . It would be interesting to see if the defect scaling exponent depends 
on the temperature range in the same way as the exponent v. In the experiment, the 
defect density would correspond to residual singlets in the final state which can be 
computed by measuring the total magnetization of the system immediately after the 
quench. We note that for these dimer systems, it would be necessary to take special 
care to achieve the criterion of long relaxation time mentioned earlier. 

1.4 Quantum communication 

In this section, we demonstrate that a properly engineered non-adiabatic dynamics 
may lead to larger fidelity and higher speed for the transfer of a qubit through a 
system. For this purpose, we begin with a Heisenberg spin- 1/2 chain described by a 
generic time-dependent Hamiltonian 

H = -Mt) E(5f5j + 5j'5p + A(t)Y^S^S) + B[t) ^5^. (1.69) 

We assume that the spin system is on a ring with sites. We start with the initial 
ground state being ferromagnetic, and denote this state by |G). At the start of the 
procedure of qubit transfer, we put a state cos(0 /2) | "f ) + sin(0 /2) exp(!0) | j) at the 
r* site of the chain. Thus the initial state of the system at the start of the evolution is 

ma 

|V/in> = cos(0/2)|G) +sin(0/2y^|r), (1.70) 

where \r) denotes the state of the spin chain with one flipped spin at the site r. We 
now consider the evolution of this state under a time-dependent Hamiltonian H. The 
specific form of the interaction need not be specified at the moment. Since the total 
spin is a conserved quantity {\^iS'f,H] = 0), the state of the system at time t becomes 

|V/(0) - cos(0/2)|G) +sin(0/2)e'^£/„,(f)|n), 

n 

where /„,(f) = {n\e-'f'"^''^'"'\r). (1.71) 
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Since the idea of communication through the chain involves performing measure- 
ment on the state at site site, we would Uke to compute the reduced density matrix 
of this site at time f . To this end, we write the wave function 

|V/(f)) =cos(0/2)|G) +sin(0/2y* ^/„,(f)|«) +sin(0/2yV.v-(f)l^), (1-72) 

where the first line of the last equation is the contribution from all terms to |i^(f)) 
where the spin in the site is t. Note that for normalization of the wave function, 
one needs 

cos2(0/2) + sin2(0/2) ^ |/„,(f)P = 1 - |/,,(f)psin2(0/2). (1.73) 

Using ( 11.721 1 and ( 11.731 1, one find that the reduced density matrix for the s* site of 
the system is 

P.v(0 - (1 - |/v,-(Ol'sin2(0/2)) I t)(T I + psin2(0/2) U)(i | 

+ ^ {e'^U{t)\ iXT I + e-'^f:,{t)\ TXi I) • (1.74) 
The fidelity of the state transfer at the given time t is thus defined as fTSl 

Fit) ^ ^JdQ{xi/Mt)\¥m) 
2 6 3 

where the integration is over the Bloch sphere involving 9 and 0. Thus to obtain 
the fidelity of a state transfer we need to obtain the matrix elements fsr{t)- To do 
this, we note that since the Hamiltonian in ( |1.69l l conserves the z component of the 
spin, an arbitrary time-dependent dynamics always restricts the system to lie within 
the subspace of one flipped spin. This allows us to write the wave function after an 
evolution through a time t to be 

\<p{t)) ^Y.^„{t)\n) ^Y.^k{t)\k), (1.76) 

n k 

where the real space basis \n) and the wave number space basis \k) are related 
by \n) = I^;tsxp(— for a chain with a periodic boundary condition. The 
Schrodinger equation for |0 (f )) now leads to the following equation for Cj^it) 

ic,{t) = {2J{t)co^{k) + i[4(f) +2B(r)]^ c,(f), (1.77) 

where we have neglected factors of 1 /N (N being the chain length which approaches 
infinity in the thermodynamic limit) in the expression for j3 {t). These equations are to 
be solved with the boundary condition c„(f = 0) = 5„r. This equation has a straight- 
forward solution 
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Q,(f) e-''(2«(0';os(<:)+j3(f))^ 
where a{t) ^ J' J{t')dt' , and I5{t) ^ J' ^[A{t') +2B{t')]dt' . (1.78) 
Using ( 11.781 1. one gets 

U{t)^(s\e-'f'"i'')^''\r)^{s\m 

= J^e'ilKr-s)+2a{t)cos{k)+P{t)]^ (2 79-) 

k 

For an infinite chain, the momentum sum can be converted to an integral and exactly 
evaluated to yield 

Mt) = Jr-si2a{t)) e-'l^^'\ (1.80) 

From this result, we note the following points. First, we need to choose a time when 
we shall perform a measurement on the state. This time, fo, is chosen so as to max- 
imize the fidelity of the state transfer. In the present model, this occurs at the time 
to when the argument 2a(fo) of the Bessel function approximately equals r~s. This 
suggests that one can reach the maximum fidelity (i.e., maximum 1/1, (0 1 ^i^^ maxi- 
mum Re[ fsr{t)] which requites a separate adjustment of the phase factor) for a given 
separation r — i at a much shorter time for a suitable non-adiabatic dynamics. This 
ensures faster communication through the channel. Note that by choosing an appro- 
priate form of 7(f), the communication can be made exponentially faster compared to 
adiabatic dynamics since we may ramp up the effective instantaneous velocity so that 
a given separation r — s is reached at a much shorter time. Second, the non-adiabatic 
dynamics gives us an additional handle on the phase and hence the real part of /s,-(f )■ 
Thus one can adjust the phase using a user-chosen classical control parameter (such 
as frequency in the case of AC dynamics) to obtain maximum fidelity for a given 
|/sr(f)|. Finally, it is straightforward to generalize the derivation of fsr{t) to higher 
dimensions. The result for a 2D system is 

/sr - 7,-.s,(2a(f)) y„,_,,,(2Ma(f)) e-'P('\ (1.81) 

where jj, is an anisotropy parameter which signifies the relative strengths of couplings 
of the Sx and Sy terms in the two orthogonal spatial directions. For jJ. — 1, i.e., the 
isotropic case, we find that the fidelity is maximized when propagation takes place 
along the diagonal. But in general, the angle of maximum propagation is a function 
of /i and this can in general also be controlled. A similar analysis can be easily 
extended to higher dimensions; however, as can be seen from (II.8II 1. the fidelity of 
the qubit transfer using this method rapidly decays with increasing dimensions. 

Thus we find, via a simple analysis of a Heisenberg spin model with a time- 
dependent Hamiltonian, that both the fidelity and the speed of quantum communi- 
cation may be improved by using suitable non-equilibrium dynamics. We have also 
shown that such a procedure can lead to direction specific state transfer in higher di- 
mensional spin systems. Since engineering such time-dependent Hamiltonians have 
become an experimental reality, this might, in principle, provide a realizable way for 
faster communication of qubits in future experiments. 
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1.5 Discussion 

To summarize, we first discussed the response of a system of interacting bosons in 
a ID optical lattice to a sudden change in a harmonic trap potential. The system 
can be mapped to a system of dipoles described by an Ising order parameter After 
the sudden shift, the order parameter oscillates in time; the amplitude of oscillations 
depends on the initial and final trap potentials. We then considered an infinite range 
ferromagnetic Ising model in a transverse magnetic field; we studied what happens 
after the field is changed suddenly. Once again, the variation of the order parameter 
(the magnetization in this problem) with time depends in an interesting way on the 
initial and final fields and the system size. 

Next, we considered what happens when a system is taken across a quantum crit- 
ical point or a critical surface in a non-adiabatic way which is governed by a quench 
time T. This leads to the production of defects; the density of defects scales as an in- 
verse power of T, where the power depends on the dimensionalities of the system and 
the critical surface and the critical exponents z and v. This was illustrated by consid- 
ering the Kitaev model which is an exactly solvable spin- 1/2 model defined on the 
honeycomb lattice; this can be solved by mapping it to a system of non-interacting 
Majorana fermions using a Jordan- Wigner transformation. We then considered the 
effect of taking a system across a QCP in a non-linear manner at time f = 0; the 
non-linearity is parametrized by an exponent a. We found that two different things 
happen depending on whether the system passes through the QCP at f = or at a 
non-zero value of f. In the former case, the power appearing in the scaling of the 
defect density with T also depends on a; in the latter case, the power is the same 
as in a linear quench (corresponding to a = 1), but the effective quench time Tgff 
depends on a. These ideas are illustrated by considering two models in ID, namely, 
a ID version of the Kitaev model, and the Ising model in a transverse magnetic field; 
both of these can solved by mapping them to systems of non-interacting fermions 
by a Jordan- Wigner transformation. We then discussed some experimental systems 
where our results for the defect scaling can be checked. 

Finally, we used some Heisenberg spin- 1/2 models in one and two dimensions 
to discuss how a qubit can be transferred across the system. In particular, we exam- 
ined how the speed and fidelity of the transfer can be maximized by choosing the 
couplings in the Hamiltonian appropriately. 

Before ending, we would like to mention two possible extensions of the work 
discussed here. First, it would be interesting to study whether the defects studied in 
Sect. I1.3l have any non-trivial topology associated with them. If they are not topo- 
logical, it would be interesting to find other ways of changing the parameters in the 
Hamiltonian in order to produce defects which do have a topological character It is 
known that topology can affect defect production in a profound way |56|. Secondly, 
the defect density discussed in Sect. 1 1.31 follows from the density matrix of a single 
site obtained by integrating out all the other sites of the system. It is interesting to 
compute the two-site density matrix and use that to obtain various measures of two- 
site entanglement. This has been studied recently both for a sudden quench ||57l and 
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a slow quench ||58l through a QCP. Other kinds of entanglement produced due by a 
quench have also been studied ||59l . 
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